package org.gis;

public class FCDTools {
	
//	public final static double PI 	= 3.14159265354;
	private final static double D2R 	= 0.017453 ;
	private final static double a2 	= 6378137.0;
	private final static double e2 	= 0.006739496742337;
	
	public static double Distance(double lat1, double lng1, double lat2, double lng2) {
		if(lng1 == lng2 && lat1 == lat2 ) {
            return 0.0;
        } else {
            double fdLambda =(lng1 - lng2) * D2R;
            double fdPhi = (lat1 - lat2) * D2R;
            double fPhimean = ((lat1 + lat2) / 2.0) * D2R;
            double fTemp = 1 - e2 * (Math.pow (Math.sin(fPhimean),2));
            double fRho = (a2 * (1 - e2)) / Math.pow (fTemp, 1.5);
            double fNu = a2 / (Math.sqrt(1 - e2 * (Math.sin(fPhimean) * Math.sin(fPhimean))));
            double fz = Math.sqrt (Math.pow(Math.sin(fdPhi / 2.0), 2) +
            		Math.cos(lat2 * D2R) * Math.cos( lat1*D2R ) * Math.pow( Math.sin(fdLambda / 2.0),2));
            fz = 2 * Math.asin(fz);
            double fAlpha = Math.cos(lat2 * D2R) * Math.sin(fdLambda) * 1 / Math.sin(fz);
            fAlpha = Math.asin (fAlpha);
            double fR = (fRho * fNu) / ((fRho * Math.pow( Math.sin(fAlpha),2)) + (fNu * Math.pow( Math.cos(fAlpha),2)));
            return fz * fR;
        }
	}


	public static void main(String[] argc) {
		System.out.println((long)Distance(116.3510240,39.9510340,116.3405300, 39.9383650));
		System.out.println((long)GetDistance(116.3510240,39.9510340,116.3405300, 39.9383650));
	}

	private static double EARTH_RADIUS = 6378137.0;
	private static double rad(double d)
	{
	    return d * Math.PI / 180.0;
	}

	public static double GetDistance(double lat1, double lng1, double lat2, double lng2)
	{
	    double radLat1 = rad(lat1);
	    double radLat2 = rad(lat2);
	    double a = radLat1 - radLat2;
	    double b = rad(lng1) - rad(lng2);
	    double s = 2 * Math.asin(Math.sqrt(Math.pow(Math.sin(a/2),2) + 
	     Math.cos(radLat1)*Math.cos(radLat2)*Math.pow(Math.sin(b/2),2)));
	    s = s * EARTH_RADIUS;
	    s = Math.round(s * 10000) / 10000;
	    return s;
	}



}